Condensation in zero-range processes on inhomogeneous networks 
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We investigate the role of inhomogeneities in zero-range processes in condensation dynamics. We 
consider the dynamics of balls hopping between nodes of a network, and find that the condensation is 
triggered by the ratio ki/k of the highest degree k\ to the average degree k. Although the condensate 
takes on the average an extensive number of balls, its occupation can oscillate in a wide range. We 
show that in systems with strong inhomogeneity, the typical melting time of the condensate grows 
exponentially with the number of balls. 
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I. INTRODUCTION 

Zero-range processes [l|, d, H, 0, H, H[ have attracted 
attention of many researchers since they provide an ex- 
actly solvable example of far-from-equilibrium dynamics 
and of condensate formation. Many questions concerning 
the dynamics of the model can be addressed and solved 
analytically. It is known that a zero-range process has 
a steady state and that this state is described by the 
partition function of the balls-in-boxes (B-in-B) model 
0, [H, also called the Backgammon model. The B-in-B 
model has two phases: a fluid and a condensed one, sep- 
arated by a critical point at which the system undergoes 
a phase transition and the condensate is formed. Unlike 
the Bosc-Einstcin condensation, the B-in-B condensation 
takes place in real space rather than in momentum space. 
Therefore it mimics such processes like mass transport 
[l[ , condensation of links in complex networks 0, Q or 
phase separation [l(| H3 • 

The zero-range process (ZRP) [m discussed in this 
paper describes a gas of identical, indistinguishable par- 
ticles hopping between the neighboring nodes of a net- 
work. The state of such a system is characterized by the 
topology of the network, which is fixed during the pro- 
cess, and by the particle distribution which is given by 
the occupation numbers of particles {to^} on all nodes 
i = 1, . . . , N of the network. The total number of parti- 
cles M = mi + m,2 + . . . + m^v is conserved during the 
process. The zero-range dynamics is characterized by 
the outflow rates it (to) of particles from network nodes, 
which depend only on the occupation number m of the 
node from which the particle hops. We shall assume that 
these ultra local hopping rates are identical for each node. 
The stream of particles outgoing from a node is equally 
distributed among all its neighbors. So if we denote by 
ki the number of neighbors of the node i, called also its 
degree, then the hopping rate from the i-th node to each 
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of its neighbors is equal to j-u(rrii). The outflow rate 
it (to) is a semi-positive function which is equal to zero 
for m = 0. For a given graph the function it(m) entirely 
defines the dynamics of the system. 

We consider here the ZRP on a network being a con- 
nected simple graph. In this case, the ZRP has a unique 
steady state, in which the probability P(mi, . . . , tojv) of 
finding the distribution of balls {toi, . . . , to^t} factorizes 
into a product of some weight functions |>i(m,) for in- 
dividual nodes, except that there is a global constraint 
reflecting the conservation of particles. 

On a A:-regular network, that is when all node degrees 
are equal to k, all weight functions Pi(m) are identical, 
Pi (to) = p(m), and thus the probability P is invariant 
under permutations of the occupation numbers. When 
the density p = M/N of balls per node exceeds a cer- 
tain critical value p c depending on the functional form 
of p(m), a single node attracts an extensive number of 
balls called the condensate. The relative occupation of 
that node does not disappear in the thermodynamic limit 
N, M — > oo, with fixed p. The larger the density p, the 
larger is the number of balls in the condensate. In other 
words, the system undergoes a phase transition at p = p c 
between the fluid (low density) and the condensed (high 
density) phase. The permutation symmetry of {to,} is 
respected in the fluid phase while it is broken in the con- 
densed phase where one node becomes evidently distinct 
from the N—l remaining ones. The symmetry of the 
partition function reduces to the subgroup of permuta- 
tions of N — 1 occupation numbers. This mechanism has 
been extensively studied in the B-in-B model 0, H, EH • 
The value of the critical density depends on the weight 
function p(m) which can be translated to the asymptotic 
properties of u[m) for to — > oo. If u(m) tends to in- 
finity then also p c — oo and the condensation does not 
occur regardless of the density of balls. The system is in 
the fluid phase for any finite density p. Intuitively, this 
means that there exists an effective repulsive force pre- 
venting a node from being occupied by many balls and 
they distribute uniformly on the whole graph. On the 
contrary, if u(m) — * 0, the critical density is p c — and 
therefore the system is in the condensed phase for any 
p > 0. The larger the number of particles on a node, the 
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smaller are the chances for balls to escape from it since 
u(m) becomes very small for large m. This can be seen as 
the existence of an effective attraction between particles. 

The most interesting case is when u(m) goes to some 
positive constant Uoc with m — * oo. One can show that 
the probability distribution P{m\, . . . , m^) does not de- 
pend on Uoo [131 ] but on how fast u(m) approaches the 
constant value. Therefore without loss of generality we 
can choose Mqo = 1 and concentrate on the asymptotic 
behavior having the form: u(m) — 1 + b/m, with b being 
some positive number. If < b < 2, then the critical 
density p c is infinite. The effective attraction between 
balls is too weak to form the condensate. However, if 
b > 2, then p c has a finite value. In this case the attrac- 
tion is strong enough to trigger the condensation above 
the critical density p c . 

So far we have discussed the criteria of the conden- 
sation on fc-regular networks, where it arises as a result 
of the spontaneous symmetry breaking. All those facts 
are well known [l^, [3] • On the other hand, the permu- 
tation symmetry can be explicitly broken if the weight 
functions Pi(m) are not identical for all nodes. For in- 
stance, this happens when the network on which the pro- 
cess takes place is inhomogeneous, that is when the de- 
grees k\ , . . . , kpi vary. A particularly important example 
are complex networks [151 ]. for which the distribution of 
degrees has usually a long tail, and thus there are many 
nodes with relatively high degree. They are, however, not 
easy for analytical studies although some predictions are 
possible [l6| . Below we shall argue that to gain some in- 
sight into the static and dynamical properties of the ZRP 
on such networks it is sufficient to study some simplified 
models. 

In the remaining part of the paper we consider only 
the most favorable situation when the weight of only 
one node is different from the remaining ones. Such 
an inhomogeneity of the weights can be introduced ei- 
ther by an inhomogeneity of the outflow rates Ui(m) or 
by an inhomogeneity of the degree distribution. We fo- 
cus here on the latter situation, when a graph has one 
node of degree k\ which differs from all remaining de- 
grees &2 = • • • = fcjv = k. As we shall see, in this case 
the quantity /i = log(fci/fc) plays the role of an external 
field breaking the permutation symmetry. 

In particular, we discuss the dynamics of the conden- 
sate on such inhomogeneous networks. This is a rela- 
tively new topic and, in contrast to the stationary prop- 
erties, less understood. Although the emergence of the 
condensate has been studied for homogeneous and in- 
homogeneous systems [1 2| | and numerically for scale-free 
networks in [lg], studies of the dynamics of an exist- 
ing condensate are rare fl7j ]. For instance, one question 
which may be asked is what is the typical life time of 
the condensate, that is how much time does it take to 
"melt" the condensate at one node and rebuild it at an- 
other node. To provide an answer to this problem is the 
main goal of the present article. 

The rest of the paper is organized as follows. In Sec. II 



we discuss the static properties of the ZRP on inhomo- 
geneous networks. We consider some particular graph 
topologies: fc-regular graphs, star graphs and fc-regular 
graphs with a single inhomogeneity introduced by a ver- 
tex of degree fci > fc. Their advantage is that all calcu- 
lations can be done exactly or at least with an excellent 
approximation. In all cases we calculate the effective oc- 
cupation number distribution 7r(m) and use it to derive 
information about the condensate dynamics. In Sec. Ill 
we derive analytic expressions for the life time of the 
condensate. We concentrate on the role of inhomogene- 
ity and typical scales at which it becomes relevant. All 
analytical results are cross-checked by Monte Carlo sim- 
ulations. The last section is devoted to a summary of our 
results. 



II. STEADY STATE - STATICS 

A zero-range process on a connected simple graph 
has a steady state with the following partition function 
Z(N,M,{ki}) O: 

M M JV 

rni— niN—0 i—1 

where &y denotes the discrete delta function and the 
weight function p(m) is related to the hop rate u(m) 
through the formula: 

m ^ 

PH=n- r r, p(0) = l. (2) 
- LJ - u[n) 

We shall denote Z(N, M, {h}) in short by Z(N, M), hav- 
ing in mind its dependence on the degrees. The partition 
function ([1]) contains the whole information about the 
static properties of the system in the steady state. The 
only trace of graph topology in the formula is through 
the nodes degrees. The dynamics, however, depends also 
on other topological characteristics but they become im- 
portant only in refined treatments. In a sense, the degree 
sequence is the first-order approximation also for the dy- 
namics. 

The probability P(mi, . . . , mjv) of a given configura- 
tion {rrii} reads: 

1 N 

P(m 1 ,...,m N ) = — - — Wpirrii)^ 
1 N 

= Z(ATM)Il^ (mi) ' (3) 

where we have defined "renormalized" weights pi(nii) — 
p(rrii)k™' i , being now node-dependent. The most impor- 
tant quantity characterizing the steady state is the proba- 
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bility 7Ti(m) that the i-th node is occupied by m particles: 



T.( m >) 



mi mi_im^i mjy 



•••X/- p ( mi '*--' miv ) x 



EN 



Zi(JV-l,M-mi)_ , . ... 

> M = z(aTm) (4) 



where Z;(A — 1,M — m) denotes the partition func- 
tion for M — m particles occupying a graph consisting 
of N — 1 nodes with degrees {fci, . . . , fcj_i, fcj+i, . . . , fcjv}- 
We shall call Pi(m) "bare" occupation probability while 
7Tj(m) "dressed" or effective occupation probability of the 
node i. We also define the average occupation probabil- 
ity: 



7r(m) = (1/N)^2n{m). 



(5) 



For a fc-regular graph, that is for fcj = fc, the occupation 
probability 7Tj(m) = 7r(m) is the same for every node 
and all the formulas above reduce to those discussed in 
Q. In general, the partition function can be calculated 
recursively: 



us as a reference point. Then we consider the particular 
example of a star topology as the simplest example of 
a single defect and finally an arbitrary fc-regular graph 
with a singular node k\ > k. The system has now the 
following weights: pi(m) — k™p(m) for the singular and 
Pi(m) = k m p(m) for the regular nodes. They differ by an 
exponential factor |>i(m)/p,(m) = (fci/fc) m = e*"™ where 
fi = log(fci/fe) > 0, which clearly favors the situation in 
which the singular node has much more particles than 
the regular ones. To make things as simple as possible, 
and to concentrate on the effect of inhomogeneity we as- 
sume that the outflow rate u(m) = 1 is constant and 
independent of m. In this case p(m) = 1 is also con- 
stant, simplifying calculations. All other functions with 
the asymptotic behavior u(m) — > 1 would lead basically 
to the same qualitative behavior. This is because in this 
case pirn) would have a power-law tail which is much less 
important for the large m-behavior than the exponential 
factor e^ m introduced by the inhomogeneity. We shall 
briefly comment on this towards the end of the paper. 



A. fc-regular graph 



Z{N,M,{h,...,k N }) = 



N-l 



m N 
M 



m 1 ,...,m J v_i 



PN(m N )Z(N — 1, M — m N , {ki, . . . , fcjv-i})- 



miv— 



(6) 



For N = 1 the partition function simply reads 
Z(l,M, ki) — pi(M). The recursive use of the formula 
([6]) allows one to compute the partition function within 
a given numerical accuracy. Using this method we were 
able to push the computation as far as to N of order 
500. For identical weights, one can find a more efficient 
recursion relation by splitting the system into two having 
similar size, which allows to study much larger systems. 
The computation of the partition function can be used 
together with Eq. ([J| to determine numerically the node 
occupation distribution 71^ (m). This gives an exact re- 
sult with a given accuracy and is more efficient than the 
corresponding Monte Carlo simulations of the ZRP. The 
dynamics, however, is not accessible in this way. 

As mentioned in the introduction we are going to ex- 
amine the effect of topological inhomogeneity on the 
properties of the ZRP. We shall consider an almost fc- 
regular graph with one node, say number one, having 
a degree bigger than the rest of nodes: fci > fc and 
fc2 = • • • = fc^v = k. The simplest realization of such 
a graph is a star or a wheel graph. In general, such a 
graph can be constructed from any fc-regular graph by 
a local modification. We proceed as follows. First we 
derive exact formulas for the particle distribution in a 
steady state on a fc-regular graph which shall later serve 



With the assumption p(m) — I, the partition function 
Z(N, M) from Eq. for the steady state of the ZRP on 
a fc-regular graph reads: 



M M 



mi— mjv— 



k M — <b dzz 
2m 



N 



-M- 



ME 



\m— 



We used an integral representation of the discrete delta 
function which allowed us to decouple the sums over 
mi, . . . , mpf for the price of having the integration over 
z. The sum over m can be done yielding 1/(1 — z) . Using 
the expansion: 



1 \ N ~ / TV 



1-z 



(8) 

and Cauchy's theorem we see that the contour integra- 
tion over z selects only the term with m = M from the 
integrand in ([7]), so we obtain: 



Z ieg (N,M) = k 



M( N + M-1 
M 



(9) 



Inserting this into Eq. ^ we find the occupation number 
distribution 



(M + N - m - 2\ . (M + N - 1 
7r(m) = m{m) = ( M _ m ]/f M 



oc 



{M + N-m- 2)! 
[M-m)\ ' 



(10) 
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FIG. 1: The "experimental" distribution 7r(m) compared to 
the theoretical prediction (|10l) for regular graphs with k — 4 
and N = 20 and for M — 20 (+) and M = 40 (x) balls. 



The distribution ir(m) is identical for all nodes and inde- 
pendent on k. It falls faster than exponentially for large 
m, therefore the condensate never appears. 

In particular one can apply these formulas to the com- 
plete graph which is just a /c-regular graph with k = 
N — 1. In Fig. Q] we see the comparison between the the- 
oretical expression (fTU)) and results of numerical Monte 
Carlo simulations for a 4-regular graph with jV = 20 
nodes and two different numbers of balls, M. The sim- 
ulations of the ZRP were organized in sweeps consisting 
on N steps each. In a single step a node was chosen at 
random and if it was non-empty a particle was picked 
and moved to a neighboring node. For each graph the 
process was initiated from a random distribution of par- 
ticles. After some thermalization, measurements of w(m) 
were done on 10 4 configurations generated every sweep. 
The results were averaged over these configurations and 
then over 5 x 10 4 independent graphs drawn at random 
from the ensemble of fc-regular graphs. 



B. Star graph 

We consider first a special case of a single inhomogene- 
ity graph, namely the star graph having N — 1 nodes of 
degree k% = . . . = fcjy = 1 connected to the central node 
with k\ = N — 1. The partition function Z stSi r(N, M) is 

oo oo 

Z staT (N,M) = J2 ■ ■ E 5 £f =1 -*,M(^-l) mi 

mi— mjv— 



m= 



M — m 



as follows from Eq. ©. It is convenient to change the 
summation index from m to j — M — m which can be 
interpreted as a deficit of particles counted relatively to 



the full occupation: 

Z st3X {N, M) = (iV - 1) M jT(JV - 1)-* ( N + J 



j=0 



(12) 



Let us assume that N ^> 1. The summands in the last 
expression are strongly suppressed when j increases so 
the sum can be approximated by changing the upper limit 
from M to oo. We obtain 



Z stai (N,M) - (N-l) M J2[i 



-1 



-TV 



-(JV-1) 



(7V-1) 



1 - 



= (Jv-i) 



M 



N - 1 
N-2 



1 

N- 1 

JV-l 



l-JV 



(13) 



Using Eq. (|4]) and the partition function Z leg calculated 
in Eq. © of the previous subsection we can determine 
the distribution of particles m at the central (singular) 
node, 



Z star (7V, M) 



(JV-1) 



m _ M i M + N — m — 2\ ( N-2 



M - m 



N - 1 



JV-1 



(14) 



Similarly, we can determine the distribution of particles 
on any external (regular) node i: 



7Tj(m) = 



(N-2)(N- l)" m 
N - 1 - (N - 1)~ M 



N-2 



(N- 1)™+!' 



(15) 



We see that 7Tj(m) decays exponentially with m while 
7Ti(m) grows exponentially for m <C M: 



7Ti(m) oc e 



l { 2M 



M (iV-1) \ 

JT ! / : \ Tj 1 M+N-2 ) _ (16) 



The growth slows down for m approaching M. At m 
M, 7Ti(m) reaches its maximal value: 



7Ti(M) = 



jV-1 



JV-1 



(17) 



which tends to e _1 when N goes to infinity. In Fig. [2] we 
compare the theoretical distributions (fT4"]) and (fTS"]) with 
Monte Carlo simulations of the ZRP for N — 20 nodes 
and M = 20, 30, 40. The agreement is very good. 

It is instructive to calculate the mean number of par- 
ticles at the central node. To simplify calculations we 
make use of the distribution = ni(M — j) of the 

deficit of particles defined above: 



5fi(j) = (l-a) i - w (-a) 



t-N, ^-/-(Jv-ir 



(18) 
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FIG. 2: The "experimental" and the theoretical (solid lines) 
particle number distributions for the star graph, for N = 20, 
and M = 20 (crosses), 30 (empty squares) and 40 (filled 
squares). The theoretical distributions were calculated ac- 
cording to the formula (|14[1 for the central node (rising curves) 
and according to Eq. (|15[) for external nodes (falling curve, the 
Monte Carlo data (points) plotted for N = 20 and M = 30). 



The parameter a — 1/(N — 1) is just the ratio of any- 
external node degree to the degree of the central node 
and measures the level of inhomogeneity. The overall 
prefactor (1 — a) 1_JV is independent of j. It is just a 
normalization constant that results from summing the 
j-dependent part of the expressions (JT8j) , 



3=0 



-(N-r 
j 



i 



(1 - a 



\N-1 ' 



(19) 



As before we changed the upper limit from M to infinity 
because a <C 1 for the star graph and hence the sum- 
mands are strongly suppressed for large j. The average 
deficit at the central node is 



(7*1 



M 
3=0 



d log 5(a) N-l 



da 



N-2' 



(20) 



as follows from Eq. (|19p. For large N it tends to one, so 
we have 



mi 



M - (ji) = M - 1. 



(21) 



We see that on average almost all balls are concentrated 
at the central node and only one ball is in the rest of the 
system. We can also determine the range of fluctuations 
around (mi) by calculating the variance. Taking again 
advantage of the generating function one finds 



(22) 



((m 1 -<m 1 )) 2 ) = ((ji-0i» 2 > 
_ d 2 logS'(e^) fN-T " 



d^ 2 



N-2 



where we used the parameter p = — log a = log (A — 1). 
For large N the result tends to one, so we can draw the 



following picture. For any JV > 1 we observe a conden- 
sation of particles at the central node regardless of their 
density p = M/N. The critical density is equal to zero 
and the system is always in the condensed phase. The 
condensate residing at the central node contains M — 1 
particles, with very small fluctuations, while the other 
nodes are almost empty. 



C. Single inhomogeneity 

A very particular property of the star graph is that 
the inhomogeneity increases with its size. Therefore, it 
is interesting to consider the situation when the inho- 
mogeneity a — k/ki is arbitrary and independent of N. 
The single inhomogeneity graph we consider here has one 
node of degree k\ and N — 1 nodes of degree k. Again, 
we assume that k± > k. The partition function (fTJ) takes 
now the form 



M 



Z inh (N,M) = J2 



k™ 1 



mi— 



M 

E 

77?2 ,. ■ .,mjv— 



(23) 

The sum over ma, . . . , tun is equal to the partition func- 
tion Z rcg (N - 1, M - mi) given by Eq. ©. The whole 
formula looks almost identical to that for the star graph 
except that now the degree k± does not need to be much 
greater than k and therefore the substitution of M by 
oo has to be done carefully in a manner incorporating 
finite-size corrections. As before, we first change vari- 
ables from mi to j = M — mi- Using Eq. (|19p we can 
cast the formula (|23|) into the following form: 



M 



Z inh (N,M) 



fcf [(1 - a 



N + j-2^ 
j 

I 1 -* - c(M)] . (24) 



The correction c(M) is equal to the sum over j from M+l 
to infinity. It corresponds to the surplus which has to be 
subtracted from the infinite sum represented by the first 
term in square brackets. In the limit M — ► oo it can be 
estimated as follows: 



c(M)= £ 



j=M+l 



N+j-2 



1 



(N-2)\J M 



dje 



Hi) 



(25) 



where 



F(j) = j log a + log ((N + j - 2)1) - logO'O- (26) 

Using Stirling's formula we can calculate the integral |25|) 
by the saddle-point method. Taking into account only 
leading terms we have 



dje 



Hi) 



,F(3.) 



2F»(j.) 



fc (M-j;)V-^'U) » (27) 
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where erfc denotes the complementary error function, 



erfc(x) = —;= 



dye' 



(28) 



and j* is determined from the saddle equation F'(j^) = 

a(N - 2) 



J* . 

1 — a 

Collecting all terms we eventually find 



(29) 



c(M) 



((JV-2)/(l-a))! Jira(N-2) 
1-a (a(N -2)/(l~ a))l 



(N — 2)! 



-erfc 



M(l - a) - a(7V - 2) 
4 ^a(iV-2) , 



(30) 



In order to keep formulas shorter we used here the nota- 
tion x\ = r(l + x) also for non-integer arguments. The 
complete partition function Z^^N, M) is given by the 
right-hand side of Eq. {23]) with c(M) given by Eq. ([30]). 
We can now calculate 7Ti(to), that is the distribution of 
balls at the singular node, 



Z Teg (N- l,M-m) 
Z inh (iV,M) 



(31) 



where Z rcg is the partition function |9]) for a regular 
graph with degree k. Using Eqs. © and (|2~4"]) we obtain 



7ri(m) = 



M + N — m — 2\ a M - ,n 

M-m ) (1 - ay~ N - c(M) ' 



(32) 



In Fig. [3] we show the theoretical ball distributions for 
graphs with k — 4, N = 20 and various M for singular 
nodes with k% = 8 and = 16, respectively, and com- 
pare them with the corresponding results obtained by 
Monte Carlo simulations. The agreement, although very 
good for the presented plots, is the better, the smaller is 
the ratio a = fc/fci. Neglecting an inessential normaliza- 
tion, we see that Eq. (|3"2"j) has the asymptotic behavior 



7T| .(III ) \ { y 



^7 + ^-771-2 
M — TO 



where 
G(to) = 



M ■ 



• exp(G(m)), 

(33) 



AT - m - - ) log(M + iV - m - 2) - 



to log a - ( M - to + - ) log(M - to). (34) 



The number of particles of the condensate can be esti- 
mated using the saddle point equation G'(to*) = for 
to* > 0. Neglecting terms of order 1/M 2 we find 



to* = M 



1 — a 



(N-2). 



(35) 




FIG. 3: The distribution of balls at the singular node for 
graphs with k = 4, N = 20, ki = 8 (top) and ki = 16 (bot- 
tom). The total number of balls is M = 20,40 and 80 from 
left to right curve, respectively. Points represent numerical 
data while solid lines show Eq. (132 p . 



Alternatively one can calculate the number of particles 
of the condensate as the mean of the distribution 7ti(to). 
Adapting the same trick as in the previous section, 



(mi) = M — ( ]1 )=M-a 



dlog5(q) 
da 



M 



1 — a 



(N - 1) M TO* 



(36) 



as follows from Eqs. (fT9"]) and ([2"0"]) . The criterion for con- 
densation is that the central node contains an extensive 
number of balls. In the limit N, M — > oo and fixed den- 
sity p = M/N it amounts to the condition (toi) > 
leading to the critical density 



Pc 



1 



(37) 



The condensation takes place when p > p c exactly like 
in the Single Defect Site model [l2j • The critical density 
decreases with decreasing ratio a — k/k\ or, equivalently, 
with increasing "external field" p — \og(ki/k). The sin- 
gular node attracts N(p — p c ) + p c balls on average as 
follows from Eq. (|3"6"|) . It is also easy to find that the 
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distribution of balls m (m) at any regular node falls ex- 
ponentially, 



7ij(m) 



fci 



m _ -urn 



(38) 



thus the condensate never appears on it. A regular node 
contains on average (rrii) = p c balls independently of the 
total density p of balls in the system as long as it exceeds 

Pc- 



III. DYNAMICS OF THE CONDENSATE 

Let us now turn to a discussion of the dynamics of the 
condensate. From the previous section we know that the 
condensate spends almost all time at the node with high- 
est degree. However, occasionally it "melts" and disap- 
pears from the singular node for a short while. We know 
that the probability of such an event is very small, so we 
expect the life-time of the condensate to be very large. 
Following the ideas of [I?) , let us imagine that we monitor 
only the number of particles at the singular node, which 
fluctuates in time. The temporal sequence of occupation 
numbers at this node performs a sort of one-dimensional 
random walk and can be viewed as a Markov chain. Us- 
ing a mean-field approximation one can derive effective 
detailed balance equations for the incoming and the out- 
going flow of particles for this node. The approximation 
is based on the assumption that the remaining part of 
the system is quickly thermalized, much faster than the 
typical time scale of the melting process on the moni- 
tored node. Thus the balance equations are written for 
the singular node and a single mean-field node having 
some typical properties. For this mean- field dynamics 
one can derive many quantities of interest. In particular 
it is convenient to calculate the average time r mn it takes 
to decrease the occupation number of the monitored node 
from m to n, or more precisely, the average first passage 
time for the Markov process initiated at m to pass n. 
This quantity was first derived in [17j for the ZRP on a 
complete graph with outflow rates u(m) — 1 + b/m. The 
formula derived there, 



E 



i 



M 



n+1 U (P)*(P) , 



(39) 



can be easily adapted to the case discussed in our paper 
by setting u(p) = 1 and using the distribution TTi(p) of 
the singular node in place of 7r(p) in the original formula. 
Equipped with the formula for r m „ we are in principle 
able to calculate the typical melting time r. What is yet 
missing is the condition for n at which the condensate can 
be considered as completely melted. We shall choose the 
simplest possible criterion and define the "typical" melt- 
ing time r as t to o, that is the time needed to completely 
empty the monitored node beginning from m equal to 
the average occupation of the node in the steady state. 



It was shown [17j that for the complete graph and 
u(m) — 1 + b/m, the melting time is approximately given 

by ' 



T OC (p - Pc ) b+1 M 



(40) 



where p c is the critical density above which the conden- 
sate is formed. The power-law increase with M can be 
attributed to the power-law fall of ir(m) ~ m~ & , charac- 
teristic for homogeneous systems with u(m) = 1 + b/m. 
The key point of our paper is that for inhomogeneous net- 
works the melting time does no longer follow a power-law 
but instead increases exponentially with M due to the oc- 
currence of the inhomogeneity which can be regarded as 
an external field p = log(fci/fc), breaking the symmetry. 

Before we do the calculations let us make a general 
remark about the dependence of T rnn on m and n. A 
quick inspection of Eq. (|3"9"1) tells us that significant con- 
tribution to the sum over p comes from terms for which 
tt(p), respectively iri(p), is small. As we know from the 
previous section, in the condensed phase 7Ti(p) is many 
orders of magnitude greater for large p than for small p. 
Therefore when m is of order M, and n is of order 1, the 
time T mn varies very slowly with m and, on the other 
hand, it is very sensitive to n. We thus put m = M for 
simplicity and concentrate on tmu- 



A. Star graph 

We assume u(p) = 1 as before. Inserting the expression 
(TT4"]) for the particle occupation distribution m (m) for the 
central node of the star into Eq. (f3"9"|) in place of 7r(m) we 
obtain 

r - V Vr.V ^- AM + N-l-2)\{M-p)\ 
kXih {M + N-p-2)\{M-l)V 

(41) 

From Sec. II B we know that the condensate contains 
m = (mi) M balls in the steady state and that fluc- 
tuations are very small. This justifies the choice m = M 
we made above. Changing the summation variables sim- 
ilarly as in the previous section we find: 



M-n-l 



T~Mn 



(N-2)l J2 



;(N-iy X 



^ (N-2 + r)\ 

^ V ; q\(N-2)\ 



(42) 



9=0 



In the second sum we can move the upper limit to infinity 
using exactly the same approximation as in Sec. II B: 



T~Mr, 



N -1 
N-2 



N-l 



M-n-l 



(N-2)\ 



r=0 



r\(N - l) r 
(N- 2 + r)V 



(43) 
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After a variable change r — > M - 
sum can be approximated as 



M-n-l 

E (M + A 

r=0 y 



(M-n-l- r)\ 



1 — r, the remaining 



-(N-iy 



(M 



n 



1) 



(M + AT - n - 



3)! ^ 

> r =0 



3-r)!' 
M + A - n - 3 



(M 



71 - 



1)(A-1) 



(44) 



such that we finally arrive at the formula 
(A-2)!(A" 



T Mr 



N- 1 



AT- 2 



l)M -n-l x 



M-n-l (M-n-l)! 
M-n-2 (M + AT -7i -3)!' 



(45) 



We see that the presence of (A— 1)~™ makes the time tm« 
indeed very sensitive to n. In Fig. [4] we see tmo compared 
to computer simulations. This complicated formula has 
a simple behavior in the limit of very large systems and 
for n = 0. In the limit of large M and for A being fixed, 
the time TMn grows exponentially with M, 



tmo ~ (A - 1) M = e^ M , 



(46) 



with /i = log(fci/fc) = log(A — 1), while for fixed density 
p = M/N and A — > oo it increases faster than exponen- 
tially, 



^pN log N 



(47) 



The approximate formulas (|4"6"]) and (|4"T)) can be alterna- 
tively obtained using a kind of Arrhenius law [Tt], EEl, 
which states that the average life time is inversely pro- 
portional to the minimal value of the occupation number 
distribution: 



T m o ~ l/7Ti(min), 



(48) 



where one thinks about the condensate's melting as of 
tunneling through the potential barrier in a potential 
V(m) = — log 7Ti (771) . In our case the potential V(m) 
grows monotonically with m going to zero, so the ball 
rather bounces from the wall at m = than tunnels 
through it, but the reasoning is the same. From Eq. IT4]) 
we have 7Ti(min) ~ (A — for fixed A and large M 

and we thus get again Eq. (|46|) , while for fixed density 
iv 1 (min) falls over-exponentially which results in Eq. (|4T|) . 

So far we have discussed the singular node. It is quite 
surprising that the formula (|39[) works also well for regu- 
lar nodes. If we blindly substitute 7Ti(fc) by Tti(k) we get 
the following expression for Ti mn : 



(A - 1)" - (A - l) m + (m - n)(A - 2)(A - 1) M 



(Af-2) 2 (A- 1) M "! 



(49) 



10° 



10' ; 



t- 10 6 - 



10= 



10" 



10 



12 

M 



13 



FIG. 4: The average lifetime tmo for the star graph with 
N — 10 calculated from Eq. (|45[l (solid line) and found in 
computer simulations (points). 



For to fixed, the transition time decreases almost linearly 
with 7i. The typical occupation of the regular node is 
much smaller than M, so we concentrate on m -C M, n — 
and JV>1. The approximate formula reads: 



A- 1 

A-2' 



(50) 



and grows very slowly in comparison to the life time of 
the condensate at the singular node. This linear growth 
can be easily understood as the minimal time needed for 
to particles to hop out from a regular node. One must 
remember that in practice we cannot observe transitions 
for large to, because the probability of having such states 
is extremely small as it stems from Eq. f| 15[) for 7Tj(m). 



B. Single inhomogeneity 

Next we consider the single inhomogeneity graph from 
Sec. II C. In the condensed phase the occupation m fluc- 
tuates quickly around (toi) and even if it is smaller than 
M we can assume that r mn f» TMn because the transition 
time TMm from M to to balls is very small in comparison 
to TMn ■ Therefore we shall concentrate again on TMn- 
From Eqs. (|32]) and ([39]) we have 



T Mr. 



M 

E 



M 



p-l ( 





M+N~l-2\ 
M-l ) 



p—n+1 l—p 



(M+N-p-2\ 
\ M-p ) 



(51) 



Changing variables we get 



T Mr, 



M 

E 



(M-p)\ 



=n+ 



M-p 



q=0 



i (M + N-p-2)\ 

,(M + A-p-g-2)! 
(M-p-q)\ ■ 



(52) 
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The sum over q can be approximated by an integral which 
can then be estimated by the saddle-point method. The 
saddle point is g* = a(N — 2) /(l — a) as in Sec. II C and 
therefore all calculations arc almost identical. In this way 
we obtain 



4^2 _ M \i-a) 1 1 2ira(N - 



2) 



Q 



(1~«) 2 



^ a(JV-2) ^, 

(53) 

where we used again the notation xl = T(l + x). The 
only dependence on p in this expression is through the 
factor aP . Thus to calculate TMn it suffices to evaluate 
the sum 




A I 



(M-p)l 



p— n+1 



(M + N-p-2)\ 



(54) 



Because every term in the sum is proportional to 1 /m (p) 
from Eq. (|3^|) , in the condensed phase the function under 
the sum has a minimum at the saddle point = m* . As 
this minimum is very deep, the effective contribution to 
the sum can be split into two terms: for small p -C tti* 
and for p 3> m*. The "small-p" part can be evaluated 
like in Eq. (|44p by pushing the upper limit to infinity and 
approximating the ratio of factorials by some number to 
power p. To calculate the "large-p" part it is sufficient to 
take the last two terms in Eq. namely for p — M and 
p = M — 1, because they decrease quickly. The complete 
formula for tmu is finally given by 



T~Mr, 



N — 2 



M 



(Mil 

q(jV-2) \| W (I -a) 2 
1 — a 



<2TTOt{N - 



Ml 



(M + N -2)1 



M + N- 2 X " 
a — I 



M 



x 1 — a 



M + N -2 
M 



, A/-1 



(a(N- 1) + 1) 
(N-l)l 



(55) 



In Fig. E] we compare this theoretical formula with t mo 
from numerical simulations. Equation (|55[) simplifies in 
the limit of large systems. When one allows M — ► oo 
while keeping N and a fixed, then the life time grows 
exponentially, 



tmo 



x M 

-) =e» M . 



(56) 



For p, — \og(N — I), that is for a star graph, it reduces 
to the formula (|46[) . In the limit of fixed density p — 
M /N > p c and for N, M — > oo: 



,N[- log(l- Q )+plog(p/Q)-(l+p) log(l+p)] 



(57) 



We see that the life time grows exponentially only if 
ki > k, that is for positive external field p — log(fci/fc). 



FIG. 5: Comparison between the "experimental" results 
(points) and the theoretical prediction (155[) (solid line) for 
tmo of a 4-regular graph with single inhomogeneity ki = 16. 
The graph size is N — 20. The inset compares tmo calculated 
from Eq. (f55)l (solid line) with that from Eq. (f56jl (dashed 
line), valid in the large-M limit. The prefactor in Eq. (1561 1 is 
chosen to match the M — > oo limit of the two formulas. 



As before, we can explore the limit when the single in- 
homogeneity graph reduces to a star graph. Inserting 
[i = - log a = log(iV-l) into Eq. (JS7J) we recover Eq. (|4T|) 
as the leading term for large N. 



IV. CONCLUSIONS 

In this paper we have studied static and dynamical 
properties of the condensation in zero-range processes on 
inhomogeneous networks. We have focused on the case 
where the network is almost a fc-regular graph except that 
it has a single node of degree fci larger than k. This type 
of network could be a rude prototype of inhomogeneities 
encountered in scale-free networks having a single hub 
with very high degree and many nodes of much smaller 
degrees. Indeed, from the point of view of the hub, the 
remaining nodes look as if they formed an almost homo- 
geneous system. We have shown that the distribution 
of balls 7Ti(ra) at the singular node has a maximum at 
to w N{p — p c ) where p c is the critical density above 
which the condensate is formed. The average occupation 
of regular nodes is equal to p c and the condensate never 
appears on them. However, the condensate is not a static 
phenomenon. It fluctuates and it sometimes melts and 
disappears from the singular node. Then the surplus of 
balls distributes uniformly on all other nodes. After a 
while the condensate reappears and its typical life time 
t grows exponentially like e^ M , where M is the number 
of balls and p — log(fci/fc) plays the role of an external 
field explicitly breaking the permutational symmetry of 
the system. This behavior is qualitatively distinct from 
that observed in homogeneous systems with a power-law 
distribution of balls, where r grows only like a power of 
M, and the symmetry is spontaneously broken. Thus 
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the transition [i = — > ^ > changes dramatically all 
properties of the system. 

In all above calculations we assumed for simplicity that 
the hop rate u(m) = 1 and thus for the homogeneous 
system there would be no condensation. However, in the 
case u(m) = 1 + b/m which produces the power-law in 
7r(m) for regular graphs, in an inhomogeneous system, 
apart from the condensation on the singular node, we 
would expect a second condensate on some regular node 
if the number of particles would be large enough to ex- 
ceed the critical density for the homogeneous sub-system. 
Thus we could expect the presence of two critical densi- 
ties pic and p2c and two condensates having completely 
distinct properties. We leave this interesting question for 
future research. 
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